/**
  CNOK project, Anyang Normal University, IMP-CAS
  \brief Calculate inclusive Rs=sigma_exp/sigma_th for -p of C12-16.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2022/08/27
  \date Last modified: 2022/08/27 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include <vector>
#include "TANucleus.h"

using std::vector;

static const int FNNPAR = 1; // 0: Horiuchi; 1: LenziRay
static const bool PAULI = false; // 0: off; 1: on
// static const double EK = 2100.; // incident energy: MeV/nucleon

int main(int argc, char const *argv[]){
  vector<TANucleus *> cv;
  TANucleus *t;

  // C12 //
  // n, sp, sn, cs, dcs, isp=true
  // cv.push_back(t = new TANucleus(12, 15.9567, 18.7207, 64.3297, 5.5536));
  cv.push_back(t = new TANucleus(12, 15.9567, 18.7207, 63.91, 6.55));
  t->SetBaseDir("config/C-p/C12");
  t->SetOrbit({"0p1/2_1/2-", "0p3/2_3/2-_1", "0p3/2_3/2-_2"});
  t->SetEx({2.1247, 0., 5.0203});

  // C13 //
  cv.push_back(t = new TANucleus(13, 17.5334, 4.9463, 39.4451, 5.9596));
  t->SetBaseDir("config/C-p/C13");
  t->SetOrbit({"0p1/2_1+", "0p1/2_0+", "0p3/2_1+", "0p3/2_2+"});
  t->SetEx({0., 2.723, 0., 0.95314});

  // C14 //
  cv.push_back(t = new TANucleus(14, 20.8310, 8.1764, 41.3081, 2.7340));
  t->SetBaseDir("config/C-p/C14");
  t->SetOrbit({"0p1/2_1/2-", "0p3/2_3/2-_1"});
  t->SetEx({3.7126, 0.});

  // C15 //
  cv.push_back(t = new TANucleus(15, 21.0795, 1.2181, 28.3508, 2.7842));
  t->SetBaseDir("config/C-p/C15");
  t->SetOrbit({"0p3/2_2-", "0p3/2_1-"});
  t->SetEx({0., 0.654});

  // C16 //
  cv.push_back(t = new TANucleus(16, 22.5531, 4.2503, 18.1660, 1.9852));
  t->SetBaseDir("config/C-p/C16");
  t->SetOrbit({"0p3/2_3/2-"});
  t->SetEx({0.});

  // results from other reactions and structural models //
  const initializer_list<double> c2sYSOX[] = { // SM with YSOX
    {0.6161, 2.7895, 0.3852}, {0.0044, 0.1219, 1.1720, 2.0037}, // C12-13
    {0.1343, 3.4401}, {1.5506, 1.2989}, {2.7216} // C14-15-16
  };
  const initializer_list<double> c2sWBP[] = { // SM with WBP
    {0.7229, 2.8594, 0.1548}, {0.0180, 0.1424, 1.0775, 1.9888}, // C12-13
    {0.1852, 3.5146}, {1.6084, 1.1651}, {2.8476} // C14-15-16
  };
  const initializer_list<double> c2sGSM[] = { // Gamow SM
    {0.2442, 3.2667, 0.3730}, {0.0001, 0.0490, 1.1157, 1.9021}, // C12-13
    {0.2245, 3.0142}, {1.8122, 1.0911}, {2.4052} // C14-15-16
  };
  const initializer_list<double> c2sNC[] = { // no core SM
    {0.6096, 2.7385, 0.0869}, {0.0046, 0.1088, 1.0840, 1.8796}, // C12-13
    {0.1141, 3.3002}, {1.5185, 1.1299}, {2.2656} // C14-15-16
  };
  const initializer_list<double> csTos[] = { // sigma_sp calculated by Tostevin
    {20.5982, 22.8958, 20.6527}, {19.6387, 18.5858, 20.9402, 20.5104}, // C12-13
    // FIXME: (C16: from GUESS, assuming that 
    {16.1085, 18.4598}, {17.1251, 16.9066}, {17.4489} // C14-15-16
  };
  const initializer_list<double> csAsia[] = { // sigma_sp calculated by my paras
    {36.8566, 29.6638, 29.6643}, {31.2968, 31.2973, 26.6435, 26.6292}, // C12-16
    {27.3394, 24.1925}, {22.1186, 22.1096}, {20.5903} // C14-15-16
  };

  for(size_t i = 0; i < cv.size(); i++){
    // cv[i]->SetC2S(c2sYSOX[i]);
    // cv[i]->SetC2S(c2sWBP[i]);
    cv[i]->SetC2S(c2sNC[i]);
    // cv[i]->SetC2S(c2sGSM[i]);
    cv[i]->SetCS(csTos[i]);
    // cv[i]->SetCS(csAsia[i]);
    // cv[i]->SetPauli(PAULI);
    // cv[i]->SetEk(EK);
    // cv[i]->SetFNNParOpt(FNNPAR);

    // if(cv[i]->GetA() == 16.)
      cv[i]->Rs();
  } // end for

  return 0;
} // end the main function
